clear all
macro drop _all
set more off

*ssc install excelcol

********************************************************************************
*				WORKING DIRECTORIES 
********************************************************************************

if ("`c(os)'"=="Windows") global ROOT "C:/Users/`c(username)'/Dropbox/PhD Maria Ulugbek Mapping indiv poverty" //Windows
if ("`c(os)'"=="MacOSX") global ROOT "/Users/`c(username)'/Dropbox/PhD Maria Ulugbek Mapping indiv poverty" //Mac

global origdata 	"$ROOT/2-data/1-raw"
global cultdata		"$ROOT/2-data/0-cultural data/Data"

global readydata  	"$ROOT/3-prog/CULTURE/replication/1_data"
global prog		 	"$ROOT/3-prog/CULTURE/replication/2_dofiles"
global shares		"$ROOT/3-prog/CULTURE/replication/3_shares"
global tables		"$ROOT/3-prog/CULTURE/replication/4_tables"
global figures		"$ROOT/3-prog/CULTURE/replication/5_figures"


********************************************************************************
* 			OPTIONS
********************************************************************************

* Data: GHA or MWI
global data "GHA"

* Model:
global model = 1
* 1 = baseline model
* 2 = model with geographic matching of cultural data
* 3 = model with age
* 4 = model with patrilineality // only for GHA

*----------------- DATA ----------------- *

use "${readydata}/${data}.dta", clear

*----------------- SETTINGS ----------------- *

* Iteration settings
set maxiter 100
global fgnlsiter "ifgnlsiterate(6)"

*----------------- MODEL VARIABLES ----------------- *

* Model 1: Baseline - patrilocality
*********************
if $model == 1 {
	global ethnovar 	 "patrilocal"
	global eta_extravars ""
}
* Model 2: Geographic matching
***********************
if $model == 2 {
	global ethnovar 	 "patrilocal_geo"
	global eta_extravars ""
}
* Model 3: Alternative specification, with age of women and men
***********************
if $model == 3 {
	global ethnovar 	 "patrilocal"
	global eta_extravars "avage_m avage_f"
}
* Model 4: Patrilineality
***********************
if $model == 4 {
	global ethnovar 	 "patrilineal"
	global eta_extravars ""
}

**********************************************

** alpha
global alpha1list 	"one nchild nwomen nmen urban" //  
global alpha2list 	"one nchild nwomen nmen urban" //	 
global alpha3list 	"one nchild nwomen nmen urban" // 
** beta
global betalist 	"one nchild nchild2 nwomen nmen urban" //  
** eta 
global eta1_nlsur 	"	c1 c2 c3 c4 nwomen nmen avage_k pboy urban $ethnovar $eta_extravars" //
global eta2_nlsur 	"c0 c1 c2 c3 c4 nwomen nmen avage_k pboy urban $ethnovar $eta_extravars" //
global eta3_nlsur 	"c0 c1 c2 c3 c4 nwomen nmen avage_k pboy urban $ethnovar $eta_extravars" //

** residual estimation
global list_f 		"nwomen urban $ethnovar $eta_extravars" //
global list_m 		"nmen urban $ethnovar $eta_extravars" //
global list_fm 		"nwomen nmen urban $ethnovar $eta_extravars" //
global list_kfm 	"nchild nwomen nmen  avage_k pboy urban $ethnovar $eta_extravars" //
global list_kf  	"nchild nwomen  avage_k pboy urban $ethnovar $eta_extravars" //
global list_km  	"nchild nmen  avage_k pboy urban $ethnovar $eta_extravars" //
global iv 			"lnd lnd2"
	
** residual form
global residlist1 "resid"
global residlist2 "resid"
global residlist3 "resid"

*----------------- SAMPLE SELECTION ----------------- *
** HH types
gen f =(nchild==0 & nmen==0 & nwomen>0) if nchild<. & nmen<. & nwomen<. 
gen m =(nchild==0 & nmen>0 & nwomen==0) if nchild<. & nmen<. & nwomen<. 
gen fm =(nchild==0 & nmen>0 & nwomen>0) if nchild<. & nmen<. & nwomen<. 
gen kfm =(nchild>0 & nmen>0 & nwomen>0) if nchild<. & nmen<. & nwomen<. 
gen kf =(nchild>0 & nmen==0 & nwomen>0) if nchild<. & nmen<. & nwomen<. 
gen km =(nchild>0 & nmen>0 & nwomen==0) if nchild<. & nmen<. & nwomen<. 
gen all=1

gen estsample=all
gen repsample=all

** Additional restrictions
if regexm("$data", "GHA")==1 { 
	sum maxage_k //maximum age of children
	gen adultchild = minage_a>r(max) & minage_a<18 if minage_a<. //HHs wit "adult" children, i.e. children below 18 but classified as adults	
	replace estsample=0 if adultchild==1 
}
if regexm("$data", "MWI")==1 {
}
	  
*----------------- DATA PREPARATION ----------------- *
	
** Age normalization
replace avage_m=avage_m/10
replace avage_f=avage_f/10
replace avage_k=avage_k/10
replace avage_a=avage_a/10

** HH types by number of children
gen c0=nchild==0 if nchild<.
gen c1=nchild==1 if nchild<.
gen c2=nchild==2 if nchild<.
gen c3=nchild==3 if nchild<.
gen c4=nchild>=4 if nchild<.
gen nchild2 = nchild^2

** Correction
replace pboy=0 if nchild==0
replace avage_f=0 if nwomen==0
replace avage_m=0 if nmen==0
replace avage_k=0 if nchild==0

** Logs
gen x=log(1+privexp) 
gen lnd = log(1+totinc)	
gen lnd2=lnd^2 

** Clothing budget shares
gen shexcl_k = cloth_c/privexp
gen shexcl_f = cloth_f/privexp
gen shexcl_m = cloth_m/privexp
gen shexcl 	 = (cloth_c + cloth_f+cloth_m)/privexp //	Total clothing budget share

** Extra vars
gen one=1
gen tiny=0.00000001 //1

** Patrilocality 
if regexm("$data", "GHA")==1 { 
	foreach x in patrilocal matrilocal neolocal {
		replace `x' = . if `x'_reg==. //to match samples
		sum `x'_reg if `x'<.
		gen `x'_geo = `x'_reg>r(mean) if `x'<. & `x'_reg<. //dummy ethno var
	}
}
if regexm("$data", "MWI")==1 { 
	foreach x in patrilocal matrilocal neolocal {
		replace `x' = . if `x'_gps==. //to match samples
		gen `x'_geo = `x'_gps if `x'<. //dummy ethno var
	}
}

*----------------- TRIMMING ----------------- *
foreach var in estsample repsample {
	preserve
		keep if `var'==1
		foreach v in f m kfm fm kf fm {
			if $model ==3 &  regexm("$data", "GHA")==1 & "`var'"=="estsample" {
				winsor2 x lnd if `v'==1 , replace cuts(0 99) trim 
				winsor2 shexcl_k shexcl_f shexcl_m if `v'==1 , replace cuts(0 98) trim 
			}
			else {
				winsor2 shexcl_k shexcl_f shexcl_m if `v'==1 , replace  cuts(0 99) trim 
			}			
		}		
		replace `var'=0 if missing(shexcl_k, shexcl_f, shexcl_m, x, lnd )==1
		rename `var' `var'_new
		keep hhid `var'_new
		tempfile `var'd
		save ``var'd'
	restore
}

merge 1:1 hhid using `estsampled', nogen
replace estsample=estsample_new
drop estsample_new

merge 1:1 hhid using `repsampled', nogen
replace repsample=repsample_new
drop repsample_new


* Drop observations with missing values in key variables and other irregularities
egen missing = rowmiss( ${alpha1list} ${alpha2list} ${alpha3list} ///
					${betalist} ${eta1_nlsur} ${eta2_nlsur} ${list_f} ${list_m} ${list_fm} ${list_kf} ${list_km} ${list_kfm}  ${iv} rexp_wb )

count if missing==0 & agemiss==0 & gendermiss==0
keep if missing==0 & agemiss==0 & gendermiss==0


** Drop only child HHs
drop if nmen==0 & nwomen==0

*----------------- IV Residuals ----------------- *

gen resid=.
quiet foreach v in f m kfm fm kf km { //
	count if `v'==1
	if r(N)>1 {
		reg x ${list_`v'} $iv if `v'==1 & estsample==1
		predict resid_`v', res
		replace resid_`v'=0 if `v'!=1
		replace resid=resid_`v' if `v'==1 
	}
}	
drop resid_*
	 
*----------------- ESTIMATION ----------------- *

macro list residlist1 residlist2 residlist3 alpha1list alpha2list alpha3list betalist eta1_nlsur eta2_nlsur 
disp `trim'

nlsur 	(shexcl_k= (nchild>0)*( {resid1: $residlist1} + (nchild>0)*exp({eta1:$eta1_nlsur}) / ((nchild>0)*exp({eta1:$eta1_nlsur}) + (nwomen>0)*exp({eta2:$eta2_nlsur}) + (nmen>0)*1 )*({alpha1: $alpha1list} + ( ({beta: $betalist}) * (x - ln(tiny + nchild) + ln(tiny + (nchild>0)*exp({eta1:$eta1_nlsur} ) / ( (nchild>0)*exp({eta1:$eta1_nlsur}) + (nwomen>0)*exp({eta2:$eta2_nlsur} ) + (nmen>0)*1 ) )   ) /*end x*/ )) /*end engel*/ ))  ///
		(shexcl_f= (nwomen>0)*( {resid2: $residlist2} + (nwomen>0)*exp({eta2:$eta2_nlsur}) / ((nchild>0)*exp({eta1:$eta1_nlsur}) + (nwomen>0)*exp({eta2:$eta2_nlsur}) + (nmen>0)*1 )*({alpha2: $alpha2list} + ( ({beta: $betalist}) * (x - ln(tiny + nwomen) + ln(tiny + (nwomen>0)*exp({eta2:$eta2_nlsur} ) / ( (nchild>0)*exp({eta1:$eta1_nlsur}) + (nwomen>0)*exp({eta2:$eta2_nlsur} ) + (nmen>0)*1 ) )   ) /*end x*/ )) /*end engel*/ ))  ///
		(shexcl_m= ( nmen>0 )*( {resid3: $residlist3} + (nmen >0) * 1	 				   / ((nchild>0)*exp({eta1:$eta1_nlsur}) + (nwomen>0)*exp({eta2:$eta2_nlsur}) + (nmen>0)*1 )*({alpha3: $alpha3list} + ( ({beta: $betalist}) * (x - ln(tiny + nmen)   + ln(tiny + (nmen>0) * 1			  			 / ( (nchild>0)*exp({eta1:$eta1_nlsur}) + (nwomen>0)*exp({eta2:$eta2_nlsur} ) + (nmen>0)*1 ) )   ) /*end x*/ )) /*end engel*/ )) ///
		if estsample==1,   ifgnls  /* trace */  eps(1e-06)  delta(1e-08)  vce(robust) noconst $fgnlsiter


* Convergence indicator
mat V=e(V)
gen converged=(diag0cnt(V)==0) // missing SE
sum converged

*-------------------PREDICTIONS ---------------- *
keep if repsample==1

** Resource shares
************************
global eta1_hat_exp "0"
global eta2_hat_exp "0"
foreach var of varlist $eta1_nlsur {
	global eta1_hat_exp "$eta1_hat_exp + _b[eta1_`var':_cons]*`var'"
}
foreach var of varlist $eta2_nlsur {
	global eta2_hat_exp "$eta2_hat_exp + _b[eta2_`var':_cons]*`var'"
}
predictnl share_k = (nchild>0)*exp( $eta1_hat_exp ) / ( (nchild>0)*exp( $eta1_hat_exp ) + (nwomen>0)*exp( $eta2_hat_exp ) + (nmen>0)*1 ) 
predictnl share_f = (nwomen>0)*exp( $eta2_hat_exp ) / ( (nchild>0)*exp( $eta1_hat_exp ) + (nwomen>0)*exp( $eta2_hat_exp ) + (nmen>0)*1 )   
predictnl share_m = (nmen>0  )*1 				    / ( (nchild>0)*exp( $eta1_hat_exp ) + (nwomen>0)*exp( $eta2_hat_exp ) + (nmen>0)*1 )  
replace share_k = . if nchild==0
replace share_f = . if nwomen==0
replace share_m = . if nmen==0
gen share_pk=share_k/nchild  //per-child
gen share_pf=share_f/nwomen //per-woman
gen share_pm=share_m/nmen //per-man

sum share_* if kfm==1 
sum share_* if fm==1

** Beta 
************************
global beta "0"
foreach var of varlist $betalist {
	global beta "$beta + _b[beta_`var':_cons]*`var'"
}
predictnl beta= $beta , se(beta_se)
gen beta_pass=( (beta-1.64*beta_se)>0 | (beta+1.64*beta_se)<0 ) //90% CI	
sum beta beta_pass if kfm==1
sum beta beta_pass if fm==1

** Marginal effects
************************
*** Patrilocality: difference in resource shares between patrilocal and matrilocal HHs
local etas_kfm "1 2"
local etas_fm	"2"
local nlcomiter "500"
	
foreach subgr in  kfm  fm { 
	foreach var in ${ethnovar} {
		foreach p in `etas_`subgr'' {	
			gen `subgr'_eta`p'_`var'_m=.
			gen `subgr'_eta`p'_`var'_mse=.
		}
	}
	
	forval p=0/1 {
		global ec`p'="0"
		global ef`p'="0"
		** Save mean values of observed characteristics of each group. 
		*** Kid
		quiet foreach var in ${eta1_nlsur} {
			sum `var' if  `subgr'==1 & ${ethnovar}==`p' 
			global ec`p' "${ec`p'} + _b[eta1_`var':_cons]*`r(mean)'"
		}
		
		*** Women
		quiet foreach var in ${eta2_nlsur} {
			sum `var' if  `subgr'==1 & ${ethnovar}==`p' 
			global ef`p' "${ef`p'} + _b[eta2_`var':_cons]*`r(mean)'"
		}
		
		* Number of children/women/men
		sum nchild if `subgr'==1 & ${ethnovar}==`p' 
		local nchild`p' = `r(mean)'
		sum nwomen if `subgr'==1 & ${ethnovar}==`p' 
		local nwomen`p' = `r(mean)'
		sum nmen if  `subgr'==1 & ${ethnovar}==`p'
		local nmen`p' = `r(mean)'
	}
	
	* Degrees of freedom
	count if `subgr'==1
	local df = r(N)-2
			
	** Calculate marginal effects
	*** Kid
	disp "`subgr' ethno: kids"

	cap matrix drop me
	cap matrix drop me_se 
	if "`subgr'"!="fm" {
		nlcom ((1/`nchild1')*exp( ${ec1} ) / ( 1 + exp( ${ec1} ) + exp( ${ef1} ))) - ((1/`nchild0')*exp( ${ec0} ) / ( 1 + exp( ${ec0} ) + exp( ${ef0} ))) , iter(`nlcomiter') df(`df')
		mat me = r(b)
		mat me_se = r(V)
		
		replace `subgr'_eta1_${ethnovar}_m = me[1,1]
		replace `subgr'_eta1_${ethnovar}_mse = me_se[1,1]^0.5
			
	}
	
	*** Women 
	disp "`subgr' ethno: women"
	
	cap matrix drop me
	cap matrix drop me_se 						
	if ("`subgr'"!="fm") nlcom ((1/`nwomen1')*exp( ${ef1} ) / ( 1 + exp( ${ec1} ) + exp( ${ef1} ))) - ((1/`nwomen0')*exp( ${ef0} ) / ( 1 + exp( ${ec0} ) + exp( ${ef0} ))) , iter(`nlcomiter')	df(`df')	
	if ("`subgr'"=="fm") nlcom ((1/`nwomen1')*exp( ${ef1} ) / ( 1 + exp( ${ef1} ))) - ((1/`nwomen0')*exp( ${ef0} ) / ( 1 + exp( ${ef0} ))) , iter(`nlcomiter')   df(`df')
	mat me = r(b)
	mat me_se = r(V)
	
	replace `subgr'_eta2_${ethnovar}_m = me[1,1]
	replace `subgr'_eta2_${ethnovar}_mse = me_se[1,1]^0.5

}	// end of loop by subgr


*** Other vars 
local nlcomiter "500"
local etas_kfm "1 2"
local etas_fm	"2"
local select_kfm "pboy urban"
local select_fm "urban"

foreach subgr in  kfm fm  { 
	
	foreach var in `select_`subgr'' {
		foreach p in `etas_`subgr'' {	
			gen `subgr'_eta`p'_`var'_m=.
			gen `subgr'_eta`p'_`var'_mse=.
		}
	}

	global ec="0"
	global ef="0"
	
	** Save mean values of observed characteristics
	*** Kid
	quiet foreach var in ${eta1_nlsur} {
		sum `var' if  nchild>0 & `subgr'==1 
		global ec "$ec + _b[eta1_`var':_cons]*`r(mean)'"
	}
	
	*** Women
	quiet foreach var in ${eta2_nlsur} {
		sum `var' if nwomen>0  & `subgr'==1
		global ef "$ef + _b[eta2_`var':_cons]*`r(mean)'"
	}
	
	* Degrees of freedom
	count if `subgr'==1 
	local df = r(N)-1
	
	sum nchild if `subgr'==1
	local nchild = `r(mean)'
	sum nwomen if `subgr'==1 
	local nwomen = `r(mean)'
	sum nmen if `subgr'==1 
	local nmen = `r(mean)'
	
	
	** Calculate marginal effects
	*** Kid
	if "`subgr'"!="fm" {
		foreach var in `select_`subgr'' { 
			
			disp "`subgr': kids: `var'"

			cap matrix drop k`var' 
			cap matrix drop k`var'_se 
			
			nlcom (1/`nchild')*( exp( $ec ) / ( 1 + exp( $ec ) + exp( $ef ) ) * ( _b[eta1_`var':_cons] - ( _b[eta1_`var':_cons] * exp( $ec ) + _b[eta2_`var':_cons] * exp( $ef ) ) / ( 1 + exp( $ec ) + exp( $ef ) ) )), iter(`nlcomiter') 	df(`df')		
			matrix k`var' = r(b)
			matrix k`var'_se = r(V)
			
			replace `subgr'_eta1_`var'_m		=	k`var'[1,1] 
			replace `subgr'_eta1_`var'_mse	=	k`var'_se[1,1]^0.5								
		
		} // end of nlcom kids
	}
	
	*** Women			
	foreach var in `select_`subgr'' {
		
		disp "`subgr': women: `var'"

		cap matrix drop f`var'
		cap matrix drop f`var'_se 
									
		if ("`subgr'"!="fm")  nlcom (1/`nwomen')*( exp( $ef ) / ( 1 + exp( $ec ) + exp( $ef ) ) * ( _b[eta2_`var':_cons] - ( _b[eta1_`var':_cons] * exp( $ec ) + _b[eta2_`var':_cons] * exp( $ef ) ) / ( 1 + exp( $ec ) + exp( $ef ) ) )), iter(`nlcomiter') df(`df') 
		if ("`subgr'"=="fm")  nlcom (1/`nwomen')*( exp( $ef ) / ( 1 + exp( $ef ) ) * ( _b[eta2_`var':_cons] - (  _b[eta2_`var':_cons] * exp( $ef ) ) / ( 1 + exp( $ef ) ) ) ), iter(`nlcomiter') df(`df')
		matrix f`var' = r(b)
		matrix f`var'_se = r(V)
		
		replace `subgr'_eta2_`var'_m=		f`var'[1,1] 
		replace `subgr'_eta2_`var'_mse	=	f`var'_se[1,1]^0.5								
	
	} // end of nlcom women

} // end of loop by HH type subgroups			



** Save results
************************			
save "${shares}/${data}_model${model}.dta", replace
	
	
